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ABSTRACT 

We obtain analytical expressions for the velocity anomaly due to the Rossiter- 
McLaughlin effect, for the case when the anomalous radial velocity is obtained by 
cross-correlation with a stellar template spectrum. In the limit of vanishing width 
of the stellar absorption lines, our result reduces to the formula derived by Ohta 
et al. (2005), which is based on the first moment of distorted stellar lines. Our 
new formula contains a term dependent on the stellar linewidth, which becomes 
important when rotational line broadening is appreciable. We generate mock 
transit spectra for four existing exoplanetary systems (HD17156, TrES-2, TrES- 
4, and HD209458) following the procedure of Winn et al. (2005), and find that 
the new formula is in better agreement with the velocity anomaly extracted from 
the mock data. Thus, our result provides a more reliable analytical description 
of the velocity anomaly due to the Rossiter-McLaughlin effect, and explains the 
previously observed dependence of the velocity anomaly on the stellar rotation 
velocity. 
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1. Introduction 

Among approximately 400 extrasolar planets discovered as of September 2009, there are 
more than 60 transiting planets. The transiting planets provide important information, such 
as radii and atmospheric signatures, which are unavailable from radial-velocity data alone. 
Another advantage of transiting exoplanets is that one can measure the angle between the 
stellar rotation axis and the orbital axis of the exoplanet projected onto the sky (convention- 
ally denoted by A) through the Rossiter-McLaughlin (RM) effect. The RM effect generates 
a radial velocity anomaly during a transit, due to the asymmetric line profiles that result 
from the partial occultation of the rotating stellar disk (Rossiter et al. 1924; McLaughlin et 
al. 1924). By carefully investigating velocity anomalies during transits, one can determine 
the trajectory of planets on the stellar disk and estimate A precisely (e.g., Queloz et al. 2000; 
Ohta et al. 2005; Winn et al. 2005). 

This angle is especially important in understanding the basic physical processes of plan- 
etary formation and subsequent orbital migration. According to recent planetary formation 
theories, gaseous planets orbiting within 0.1 AU of the parent star (hot-Jupiters) are sup- 
posed to have formed a few AU away from the star, and subsequently to have migrated 
inward. While standard migration mechanisms keep A ~ 0°, scenarios such as planet-planet 
scattering may change A significantly (Lin et al. 1996; Chatterjee et al. 2008; Nagasawa et 
al. 2008). The Kozai mechanism may also lead to a highly inclined orbit (Wu et al. 2007). 
Thus, the observed distribution of the angles provides an observational clue to distinguish 
between different planet formation and migration theories. 

Ohta et al. (2005) derived an analytic formula (hereafter, the OTS formula) for velocity 
anomalies due to the RM effect, by computing the first moment (intensity- weighted mean 
wavelength) of distorted absorption lines perturbatively. (Previously, Kopal et al. (1942) 
and others used the first-moment method for eclipsing binary stars to estimate the velocity 
anomalies.) The OTS formula proved to be useful in understanding the parameter depen- 
dence of the velocity anomaly and also in forecasting the error budget of the parameter 
estimate, in particular of the spin-orbit misalignment angle A. Indeed their work inspired 
Winn et al. (2005) to revisit the estimate of A precisely for the first discovered transiting 
planetary system HD 209458. For that purpose, Winn et al. (2005) generated in-transit 
stellar spectra of HD 209458, put them into the Keck analysis routine of radial velocities, 
and found that the OTS formula systematically underpredicts the velocity anomaly by 10%. 
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Intriguingly, however, Johnson et al. (2008) and Winn et al. (2008), found the OTS formula 
to provide an adequate description of the simulated results for the cases of HAT-P-1 and 
TrES-2. 

We noticed that the systematic difference between the OTS formula and the simulation 
may be sensitive to the stellar spin rotation velocity v sin (i* denotes the inclination angle 
of the stellar spin axis), which is 4.70 it 0.16 km s~^, 3.75 ±0.58 km s~^, and 1.0 ±0.6 km s"^, 
for the HD 209458, HAT-P-1, and TrES-2 systems, respectively. In addition, the actual radial 
velocity measurement algorithm does not adopt the moment method, strictly speaking. For 
instance, the analysis routines used by the HARPS (the High Accuracy Radial velocity Planet 
Searcher) and SOPHIE (Spectrographe pour I'Observation des Phenomenes des Interieurs 
stellaires et des Exoplanetes) teams involve the cross-correlation of the observed spectrum 
with a template spectrum, and the velocity anomaly is estimated from the peak of the 
cross-correlation function. The analysis routine which uses the iodine cell technique, with 
the Subaru HDS (the High Dispersion Spectrograph) or Keck HIRES (the High Resolution 
Echelle Spectrometer), involve a fit to the observed spectra, based on a template spectrum 
and the known spectrum of the iodine cell. 

For these reasons, in this paper we revisit the analytic approach to the RM effect using 
the cross-correlation method. Our analytic approach to this problem is complementary to 
the recent numerical approach by Triaud et al. (2009), and provides an analytic framework 
for understanding their results. Even though the cross-correlation method is not directly 
applicable to the analysis of the iodine cell technique, we find that our result seems to capture 
the main qualitative features of the numerical results based on that technique; specifically, 
our formula includes a term dependent on the stellar spin velocity that does not show up in 
the moment method but was empirically found by Winn et al. (2005). 

The rest of the paper is organized as follows. Section 2 describes our analytical modeling 
of stellar absorption lines distorted by the occultation due to a transiting planet (§2.1) and 
derives the RM velocity anomaly on the basis of the moment method (§2.2) and the cross- 
correlation method (§2.3). In order to check the reliability of the new analytic formula, we 
generate simulated spectra during transits and put them into the Subaru analysis routine 
of radial velocities in §3 as Winn et al. (2005) did for the Keck routine. Specifically we 
consider HD 17156 (f sini^ = 4.2 km sec~^), TrES-2 (f sini^ = 1.0 km sec"^), TrES-4 
{vsmii, = 8.5 km sec~^), and HD209458 (wsini^ = 4.5 km sec~^) systems, and find that our 
new formula reproduces the simulated data better than the OTS formula. Finally, section 4 
is devoted to a summary and further discussion of the present paper. 



-4- 



2. Analytic Model for the Velocity Anomaly 



In this section, we analytically derive the velocity anomaly due to the RM effect based 
on the two different methods. The final result for the moment method is equation (16). For 
the cross-correlation method, the key result is equation (37) when the stellar line profile and 
the rotational broadening kernel are approximated by equations (20) and (22), respectively. 
Readers who are not interested in the mathematical details of the derivations may simply 
skip this section and directly move to §3. 



The stellar absorption line profile distorted by the RM effect is schematically illustrated 
in Figure 1. During a transit, a portion on the stellar disk is occulted by a planet, and the 
flux contribution from the portion occulted by the planet, jFpianct(A — AA) is subtracted from 
the stellar absorption line, J?star(A), yielding a distorted absorption line profile in transit, 
-^transit (A). The Center of J-pianet is Doppler-shifted by AA relative to the central wavelength 
Aq of J^sta.v since the portion occulted by the planet has a stellar rotation velocity: 



along the line of sight of an observer. Because of the distortion, the absorption line profile 
becomes asymmetric, producing an apparent velocity anomaly when the radial velocity is 
estimated using one of the standard techniques. 

Let us first write a stellar absorption line profile as a convolution of an intrinsic line 
profile S{\) and a stellar rotation kernel -R(A; A^) as 



where the symbol * indicates a convolution, and the continuum level is normalized to unity 
(Fig. 2). In this section, we implicitly focus on a particular single line of a central wavelength 
Aq, and locate the center of S{\) at Aq. Then we normalize S{\) and -R(A; A^,) so that 



2.1. Velocity Anomaly During Transit 



AA 



(1) 




-^star(A) = l-5(A)*i?(A;A,.), 



(2) 




(3) 



A Voigt profile (convolution of a Gaussian and a Lorentzian; see eq. [19] below) is often 
used to approximate 5'(A), but a Gaussian profile also provides a reasonable approximation. 
We will mainly be concerned with the Gaussian approximation in § 2.3. 
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Fig. 1. — Schematic absorption line profiles of a transiting planetary system: J?star(A) is a 
symmetric line profile with respect to the central wavelength of each line Aq over the entire 
stellar disk outside the transit, while JFpianct(A — AA) is a symmetric line profile for the portion 
occulted by a planet that is shifted by a wavelength of AA with respect to Aq. The resulting 
distorted stellar line profile in transit is given by J^transit(A) = ^star(A) — ^pianct(A — AA). 

The function -R(A; A^) indicates a rotational broadening kernel and broadens the intrinsic 
stellar line at Aq over the range of Aq — A^ and Aq + A/,, where Xl/\q = v sini^/c. Its expression 
for a rigid rotating stellar disk can be computed once the limb darkening law is specified (see, 
e.g., Gray 2005). In this paper, we consider the quadratic limb darkening law in which the 
intensity of the stellar disk at a position (x, y) relative to the center of the star is expressed 
as 



In the above expression, 6 is the angle between the line-of-sight and the normal vector to the 
local stellar surface, Rs is the stellar radius, the coefficients ui and U2 are the limb- darkening 
parameters, and Ic indicates the intensity at the center of the stellar disk. 

Then we obtain 





(5) 
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Fig. 2. — Stellar absorption line profile Tst&v, which is a convolution of the intrinsic stellar line 
profile (including the natural broadening), ^(A) and the stellar rotational kernel -R(A; A^). 

where the coefficients ci, C2, and C3 are given as 

2(1 - ui - U2) 
7rAz,(l - Ml/3 - M2/6) 
^ Ui + 2^2 . . 

' 2Ai(l-ni/3-n2/6)' 

C3 = 



37rAL(l - Ml/3 - M2/6) 



Adopting equation (2) as a line profile, the line distortion due to the RM effect is 
described as follows. In the transiting planetary system, the portion occulted by the planet 
is sufficiently small. Thus the effect of the rotational broadening within the portion is safely 
negligible, and one can write 

J-pianct(A - AA) = / {1 - 5(A - AA)} , (9) 

where / is the ratio of the flux from the occulted part of the stellar disk to the total flux. 

In general, A A = VpXo/c is a function of the position of the planet and the contribution 
of wavelength shift is determined by integrating over the entire occulted portion of stellar 
disk. The precise definition of Vp in terms of the stellar rotation velocity and the position of 
planet is given by equation (A7), which is simplified as equation (A8) for rigid rotation. 
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Then, the line profile during transit is given as 

^ transit(A) = T 3tar(A) — T planct('^ ~ AA) 

= {l-5(A)*i?(A;AL)}-/{l-5(A- AA)}. (10) 

Because the continuum level 1 — / in equation (10) does not change the evaluation of the 
velocity anomaly, in what follows we simply use the expression 

^transit(A) = -5(A) * R{\- A^) + fS{\ - AA) (11) 

without losing generality. 



2.2. Estimate Based on the Moment Method 

Let us consider first the moment method, in which we estimate the wavelength shift of 
the center of a distorted line during a transit as (e.g. Ohta et al. 2005): 



transit ('^)'^'^ 

5rm = — • (12) 

^ transit ( A) (i A 



Substituting equation (10) into J?-transit(A) in equation (12), the denominator is 
-^transit(A)rfA = - / rfA'/?(A'; A^) / d\Si\-\') + f dXSiX-AX) 

-oo J ~oo J —oo J —oo 

= -! + /• (13) 
The numerator of equation (12) reduces to 

/oo POO POO POO 

A^transit(A)rfA = - / dA'i?(A'; A^) / XSiX-X')dX + f A5(A - AA)dA 
CO J — oo J — oo J — oo 

/oo POO 
R{X'; XL)dX' / (A - A')5(A - A')ciA 
oo J —oo 

POO 

X'R{X'; XL)dX' / S{X - X')dX 



oo 
oo 



+/ / (A- AA)5(A- AA)dA + / / AA5(A - AA)rfA 



oo 
oo 



/ / AA5(A-AA)rfA 

J —oo 

/AA, (14) 
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and therefore we obtain 



Srm = -J^^^- (15) 



Using equation (1), one can rewrite equation (15) in terms of the velocity anomaly as 

Av = C-— = -Vp. (16) 

^0 -L — / 

When we exphcitly write the flux ratio / and the subplanet velocity Vp clS cL function 
of the planet position, we reproduce the OTS formula (see Appendix A and B below). In 
particular, the above expression coincides with equation (25) of OTS which is the formula 
for complete transiting case without limb- darkening, if we simply replace the quantities of / 
and Vp with 7^ and fi^xsinz^. 

Note that equations (13) to (15) are derived from the normalization condition (3) of 
S{\) and R{\; Xl) and the fact that they are even functions of A, and do not depend on their 
speciflc functional forms. This implies that the velocity anomaly derived from the moment 
method depends on neither the proflle nor the width of the absorption line. Thus the results 
are independent of the stellar rotation velocity. 



2.3. Estimate Based on the Cross-Correlation Method 

2.3.1. Formulation for the Voigt profile 

Consider next the estimate of the velocity anomaly based on the cross-correlation 
method. In this case, the wavelength shift of central line Sbm is obtained by maximizing the 
cross-correlation function: 

dC{x 



dx 



= 0, (17) 

a;=ORM 



/oo 
-^star(A - x)J'transit(A)(iA. (18) 
-oo 

To proceed further, we have to specify the functional form of 5'(A) and -R(A, Xl). 
Here we adopt the Voigt proflle for the stellar line: 

S{X) = Gs(A;/5)*L(A;7), (19) 
Gs{X;P) = ^e-'^/f^', (20) 

^(^'^-^^ = lwTV2- (21) 

vr A"' + 7^ 
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In the above expressions, /3 in the Gaussian kernel characterizes the thermal and microtur- 
bulent broadening, while 7^ in the Lorentzian kernel represents the line broadening due to 
natural broadening and the Stark effect. 

While we derive equation (5) for the rotational kernel R{X; Xl) for the quadratic limb- 
darkening law, the expression is not analytically tractable. Thus we approximate it by a 
Gaussian: 

1 



R{X;XL)^Gn{X;a) 



aWTT 



(22) 



where the relation between a and Al is obtained by least-square fitting of Gaussian to 
equation (5). Their explicit relation is given in Appendix F. 



Under this Gaussian approximation for the rotation kernel, we have 

S{X) * R{X; Xl) ^ Gs(A; /5) * L(A; 7*) * Gr(A; a) = G(A; /5.) * L(A; 7.), 

where 

This significantly simplifies the analytic computation, and we now have 

^ transit (A) = " 



^ star (A) 



(23) 



(24) 



G'(A';/3jL(A-A';7.)^iA' 

/oo 
G(A';/?p)L(A- AA-A';7p)t^A' 
-00 

G'(A';/5,)L(A-A';7,)rfA', 



(25) 
(26) 



where /3p and 7^ indicate the narrow line originated from the portion of the planet, and thus 
we simply set Pp = p. 

Substituting equations (25) and (26) into equations (17) and (18), we obtain the follow- 
ing algebraic equation for 5rm: 



53/2 



a — tcje 



iac/2b 



1 + erfi 



a — ic 



53/2 



/ e-(^'-^')/^^Re|(A-zC)e^-^^/2B 



1 + erfi 



2Vb 
A-iC 



2VB 



where 



a = 2n5^M, b = 27i^p^, c = inj^ 
A = 27r(5RM -AX), B = tt^P^ + /J^), C = 27r(7, + 7^). 



(27) 



(28) 
(29) 
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and erfi(z) denotes the imaginary error function defined in terms of the error function erf(x): 



erfi(^) = -i ei{{iz) = / e'^^dt. (30) 



The derivation of equation (27) is summarized in Appendix C. 



2.3.2. Analytic formula for the Gaussian profile 

Since the Voigt profile is widely used for stellar absorption lines, equation (27) is fairly 
general. While it cannot be solved exactly for 5rm, we obtain an approximate perturbation 
solution for 5rm ^ AA. In particular, we present the result for the Gaussian profile (7^ = 
7p = 0) in this subsection, which will be compared with simulated data analysis in §3. The 
other solutions for the Lorentzian profile (/5^ = 0) and the Gaussian profile with a small 
Lorentzian wing (7^ <^ (3i,) are presented in Appendix D and E, respectively. 

In the Gaussian profile assumed here, the parameters c and C in equations (28) and 
(29) vanish. Then equation (27) reduces to 

^,-./4. ^ /_Le-^/^^ (31) 



J,3/2 ■> 



or equivalently 



. _ , ' P^ + Pp \ ^''^ '^RM ^^^^ / , ((5rm - AA)^ . 



V 



Note that / is approximately the square of the planet-to-star radius ratio, and is of 
order 10~^ for giant planets. Thus, equation (15) implies that |5rm/^A| ~ / <^ 1, and we 
can expand the exponent in equation (32) up to the linear order of (5rm: 



'5rm , ('5rm - AA)^ _ AA^ / JmA\ _ ( , J 



In general we consider the case where AA is well within the width of the stellar broadening, 
and therefore k ^ 1. Therefore we expand equation (32) up to the linear order in 5rm/AA 
as 



ix'^ + M. (34) 
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and finally obtain 



p 

As before, the above equation can be rewritten as 



- ( /-J 1 - -TT^r^vl } (36) 



or 



A.^-i ?(^!±^ V" fvll- ^ (37) 



where we replace /3p with (line width without rotational broadening), and with + 
in which a indicates the stellar rotation width. 

Equation (37) is the key result of the present paper which describes the RM velocity 
anomaly under the Gaussian approximation for the intrinsic line profile and the stellar ro- 
tation kernel. In marked contrast to equation (16) based on the moment method, equation 
(37) indicates that the radial velocity anomaly from the cross-correlation method depends 
both on the width of the line profile {(5) and on the rotation velocity (a). It recovers the 
OTS formula only when a = and the second term in the parenthesis is negligibly small. 

Indeed the presence of the cubic term in Vp is consistent with the empirical finding from 
the simulated data (Winn et al. 2005). Note that a quintic (5th order term) was needed to 
match the simulations for the rapidly rotating star HAT-P-2 (Winn et al. 2007), which would 
emerge naturally in our formulation by expanding the next-order term in k. We should note 
here that in many cases the estimate of the spin-orbit misalignment angle A is not very 
sensitive to the choice of the moment method or the cross-correlation method discussed in 
this section; when A is small and the transit impact parameter is not too close to zero, then 
A is mainly determined by the epoch of the photometric central transit when Av, and thus 
Vp, vanish. 



3. Comparison with Simulated Results Based on the Subaru Analysis Routine 

So far, we have focused on the RM velocity anomaly estimated from single line profiles 
characterized by the Voigt function. The practical data analysis, however, is much more 
complicated; the wavelength shift AA is estimated statistically from many different lines, and 
each line is not always approximated by a single Voigt profile. In addition, Subaru and Keck 
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analysis routines are not related in a straightforward way to the cross- correlation algorithm. 
This is why we create mock spectrum data during transit, compute radial velocity anomalies 
using the Subaru analysis routine, and compare the output results with our analytic formulae, 
in particular equation (37). 



In order to check the validity of our analytic formula, we create mock spectra in transit, 
and reduce the data with the analysis routine (Sato et al. 2002) for the Subaru HDS radial 
velocity measurement with I2 cell. The routine has been extensively used for the previous 
observation of the RM effect (e.g. Narita et al. 2007, 2008). 

The routine finds the best-fit value of the wavelength shift SX by modeling an observed 
spectrum Jobs (A) with the given transmission function of the I2 cell, A{X), as 



where k is the normalization factor, ^(A — 6X) is the intrinsic stellar spectrum shifted by 
6X, and IP denotes the instrumental profile. Sato et al. (2002) use the Lick-Hamilton I2 cell 
spectrum for A{X) with high spectral resolution {R ~ 400,000). The peculiar velocity and 
the Keplerian motion of the star as well as the RM velocity anomaly contribute to SX, and 
thus the former two have to be subtracted in order to extract 5rm alone. 

Once A{X) and S{X) are given, the routine simultaneously finds the best-fit values 
for k, 6X, parameters characterizing IP (usually expressed as a linear combination of several 
Gaussians), and calibration parameters of wavelength by least-square fit to Jobs (A). In reality, 
the method of Sato et al. (2002) attempts to elaborate the intrinsic stellar template spectrum 
S{X) as well from the observed spectrum /obs(A). In practice, it proceeds as follows: 

1. Choose an initial stellar spectrum S^{X) from an observational dataset or from a the- 
oretical model. 

2. Find the best-fit parameters in the right-hand-side of equation (38), and denote the 
resultant best-fit spectrum by I^{X). 

3. Compute the residual: 



3.1. Subaru Analysis Routine 



Jobs(A) = k[A{X)S{X - 6X)] * IP 



(38) 



6S{X) 



-^obs(A) 



(A) 



(39) 



and use 



S\X) 



S''{X)+5S{X), 



(40) 
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revised template for the stellar spectrum. 
4. Iterate the above steps until the given accuracy is achieved. 

In general, the above method finds the best-fit stellar template spectrum separately for 
each Jobs (A) observed at different epochs. In order to avoid the contamination caused by the 
I2 absorption lines in using only one spectrum, the set of stellar templates at different times 
are averaged and the resulting stellar spectrum S{X) is used to determine 6\ for each /obs(-^)- 

The Keck HIRES analysis routine by Butler et al. (1996) uses the stellar template 
spectrum directly obtained from the observed spectrum (without the I2 cell) by deconvolving 
the instrumental profile that is estimated based on observations of a rapidly-rotating B star 
using the I2 cell. Except for the treatment of the stellar spectrum, both routines are very 
similar, and our simulation result below is applicable to the Keck data analysis. 

3.2. Simulated Spectra 

As mock observational spectra for the Subaru analysis routine (Jobs (A) in eq. [38]), we 
create many simulated spectra during transits following the procedure of Winn et al. (2005) 
as follows: 

(i) Broaden the NSO solar spectrum (Kurucz et al. 1984) so as to the include the effect 
of stellar rotation^. We apply two different broadening kernels as discussed below. 

(ii) Compute the spectra during transit from equation (10) using the broadened spectrum 
described above for J^sta.rW, and the original (unbroadened) spectrum for J?-pianot(A — 
AA). 

(iii) Multiply J^transit(A) by A{\), the transmission function of the I2 cell (eq. [38]). 

(iv) Convolve J?-transit(A)A(A) with the instrumental profile, IP, which is approximated by 
a single Gaussian that best-fits the IP of Subaru/HDS. We call the resultant spectra 
"the mock data". 

(v) Feed the mock spectra into the Subaru analysis routine (§3.1), and compute SX. 



^Although the NSO solar spectrum is a rotationally broadened spectrum whose rotational velocity 
wsini* ~ 1.85 km sec~^, we additionally broaden the spectrum by Gaussian broadening kernels. 
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For the broadening kernel for the stellar rotation in step (i), we first adopt a Gaussian 
(§3.3.1) so as to compare the analytic formula under the Gaussian approximation, and then 
use equation (5) for more realistic comparison. 

Once the broadening kernel is specified, the remaining parameters are the flux ratio / 
and the subplanet velocity Vp = cAA/A. We repeat the above procedure at grids on {f,Vp) 
plane. We select / = 0.00, 0.004, 0.008, 0.010, 0.012, 0.016, 0.020, and assign to Vp typically 
20-30 different values ranging over several km sec~^ for each / in equal interval. Then we 
compute the velocity anomaly Av on those grids of (/, Vp) . 

Finally we apply the methodology to four existing planetary systems (HD17156, TrES- 
2, TrES-4, and HD209458), and fit the result to the following form inspired by the Gaussian 
approximation (eq. [37]): 

Av = -fvpip-qvl), (41) 

where p and q are constants fitted for each planetary system, and depend on the stellar 
absorption line profile and rotation broadening among others. 



3.3. Results 



3.3.1. Ideal case with Gaussian Rotational Kernel 



Consider the NSO spectrum additionally broadened with the Gaussian kernel (eq. [22]). 
We adopt a = 0, 2.7, and 6.8 km sec~^, ^ create mock spectra during transits for grids on 
(/, Vp) plane, and determine the wavelength shift AA. The velocity anomaly Av against Vp 
is plotted in Figure 3 for / = 0.01. Since Av is simply proportional to /, this plot can be 
easily scaled for an arbitrary value of / as long as / ^ 1. 

Different absorption lines in the NSO spectrum have different line profiles and are not 
exactly Gaussian. If the Gaussian approximation (eq. [37]) were exact, the coefficients in 
equation (41) for a particular line at wavelength Aq should be given as 



P 



1 + 



2/52 + C72 



3/2 



1 + 



2/32 



3/2 



^0 



c2(2/32 + a2) 



(42) 



Apparently they should be sensitive to each absorption line profile (its central wavelength 
Ao and width /3), and it is not clear that the real data analysis based on many different lines 



^The dispersion parameter a (or (3) has a unit of length, but we here express them in terms of radial 
velocities in order to make it easier to compare the dispersion parameters with the stellar rotational velocities. 



- 15 - 




Fig. 3. — The RM velocity anomaly against Vp for / = 0.01. symbols: results extracted from 
mock data using the Subaru analysis routine with NSO spectrum for S{X) and the Gaussian 
broadening kernel (eq. [22]). Crosses, circles, and triangles correspond to a = 0, 2.7, and 6.8 
km sec~^, respectively. The number of data points shown is suppressed to make it easier to 
see. solid lines: polynomial fit (eq. [41]) to the mock results. The OTS formula is plotted in 
dashed line for comparison (almost indistinguishable from the a = line). 



overall reproduces the polynomial form like equation (41). Nevertheless Figure 3 indicates 
that the simulated results are indeed fitted very well by equation (41). As expected, the 
OTS formula derived from the moment method agrees with the simulated result if the stellar 
rotational broadening is negligible (a = 0). 

Given the above encouraging result, we can use equation (41) to compute the velocity 
anomaly curve. In order to do so, we have to compute / and Vp as a function of time, which 
can be done using the formulae in Appendix A and B (see also OTS). Figure 4 plots the 
prediction for HD209458 in symbols on the basis of the result plotted in Figure 3. Also 
plotted for comparison is the analytic formula (37) for the Gaussian kernel adopting (3 = 2.6 
km sec~^ as the best-fit of the original solar spectrum. We here adopt Aq = 5500 A as a 
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time [h] 



Fig. 4. — The RM velocity anomaly curve for HD209458 but with varying a for the Gaussian 
broadening kernel. The central transit time is set to the origin of the time. The spin-orbit 
misalignment angle A is assumed to be 0°. symbols: predictions based on the polynomial fit 
(eq. [41]) plotted in Figure 3. Crosses, circles, and triangles correspond to a = 0, 2.7, and 6.8 
km sec~^, respectively, solid lines: analytic formula (eq. [37]) with the Gaussian broadening 
kernel based on the cross-correlation method. The OTS formula is plotted in dashed line for 
comparison. 

typical wavelength of the radial velocity analysis^. Since the actual stellar rotation speed of 
HD209458 is v sinz* = 4.70 ± 0.16 km sec~^, a ~ 4 km sec~^ is appropriate, but we consider 
(7 = 0, 2.7, and 6.8 km sec~^ so as to examine the dependence on a. 

Again we stress here that the lines in the NSO spectrum are not perfectly Gaussian, and 
the Subaru analysis routine adopts multi-parameters fitting with the least-squares method, 
instead of the cross-correlation method we assumed in the analytic treatment. Nevertheless 
the simulated result and the analytic formula (37) show good overall agreement, indicating 
the practical validity of our analytic approach. 



•^Many absorption lines of the iodine cell are located between 5000 A and 6000 A. 
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3.3.2. Application to actual exoplanetary systems 



While we find that the basic features of the simulated are explained by the Gaussian 
approximation both of the stellar line profile and the stellar rotation kernel, it still remains 
to see the extent to which a more realistic rotation kernel affects the RM velocity anomaly. 
We examine this point by repeating the simulated analysis for HD 17156, TrES-2, TrES-4, 
and HD209458 with the rotation kernel (5). 

HD17156 This transiting planetary system was discovered by Fischer et al. (2007). We 
assume the linear limb darkening law with ui = 0.6 and U2 = and the rotational 
velocity of w sini^ = 4.2 km sec~^(Narita et al. 2009).^ We perform the same simulation 
as before but now using the rotation kernel (5), and obtain 



The result is plotted in crosses in Figure 5 (a). The misalignment angle A is assumed 
to be 10°. We adopt io = 87.2° as the orbital inclination of HD17156b. 

While the rotation kernel is not Gaussian, it is interesting to see the difference between 
our analytical formula and the simulated result. For this purpose, we fit equation (5) 
by a single Gaussian and find the best-fit of a ~ Al/1.31 (see Appendix F for detail). 
The solid line in Figure 5 (a) indicates the corresponding result from equation (37) for 
f3 = 2.6 km sec~^ and a = i;sinz*/1.31 ~ 3.2 km sec~^. 

TrES-2 The central star in TrES-2 system has a small projected rotational velocity. Winn et 
al. (2008) investigated the RM effect of this system, and concluded the orbit of TrES- 
2b is prograde. They derived the empirical formula for the velocity anomaly using 
templates based on observations of similar stars. Here, we use a theoretical GO type 
template spectrum without rotational broadening (Coelho et al. 2005), instead of the 
NSO spectrum, to create the mock spectra. We then broaden the theoretical spectrum 
assuming vsini^ = 1.0 km sec~^ and the limb-darkening parameters Ui = 0.40 and 
U2 = 0.30. After inputting the mock spectra into the analysis routine, we obtain the 
following empirical formula: 



''While Fischer et al. (2007) estimated the stellar rotational velocity to be usini^, = 2.6 km sec~^, Narita 
et al. (2009) obtained the stellar rotational velocity ?;sini* = 4.2 km sec~^ by the RM analysis. We adopt 
the latter value here. 




(43) 




(44) 
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(b) 



Gaussian approximation 
equation (37) [a = 0.77 l^m sec^] - 

OTS\ 



RM velocity curve 
TrES-2 
misalignment angle A = ° 



-2 -1.5 -1 -0.5 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0.5 1 1.5 2 

time [h] time [In] 




-2 -1.5 -1 -0.5 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0.5 1 1.5 2 



time [h] time [h] 

Fig. 5. — The RM velocity anomaly curves for HD17156 (a), TrES-2 (b), TrES-4 (c), and 
HD209458 (d). Crosses indicate the simulated result using the polynomial fit (eqs.[43], [44], 
[45], and [46]). Solid and dashed lines correspond to the Gaussian formula (eq. [37]) and the 
OTS formula, respectively. The spin-orbit misalignment angle A is assumed to be 10° for 
HD17156 and 0° for other systems. We adopt 87.2°, 83.6°, 82.6°, and 86.7° for the orbital 
inclination of HD17156, TrES-2, TrES-4, and HD209458, respectively. 
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Again, in Figure 5 (b), crosses indicate the simulated result (eq. [44]). We also show 
the OTS formula and the Gaussian formula (37), in which a = t'sinz^/l.SO = 0.766 
km sec~^. We here adopt io = 83.6° and A = 0°. 

The three different lines are almost indistinguishable due to the small RM amplitude, 
but we can safely say that the OTS formula as well as the Gaussian formula (37) well 
describe the simulated result. When the stellar spin velocity is small enough, the OTS 
formula provides a good approximation to the velocity anomaly. 

TrES-4 TrES-4 is a transiting planetary system discovered by Mandushev et al. (2007). 
The parent star is a F-type star with an effective temperature of 6200 ± 75 K and has 
a relatively high rotation velocity of wsini^ = 8.5 ±0.5 km sec~^ (Sozzetti et al. 2009). 
Since the spectral type of the star is different from that of the Sun, we use theoretically 
synthesized F-type spectrum. We then broaden the theoretical spectrum assuming a 
quadratic limb-darkening law (eq. [4]) with ui = 0.46 and U2 = 0.31 (Claret 2004), 
and obtain 



Figure 5 (c) shows the result. The misalignment angle A and the orbital inclination io 
are assumed to be 0° and 82.6°, respectively. The solid line in Figure 5 indicates the 
Gaussian formula with a = t;sinz^/1.32 ^ 6.4 km sec~^. 

HD209458 Finally, we compare the Gaussian formula with the simulated radial velocity 
anomaly of HD209458 system. While in Figure 4 the Gaussian broadening kernel is 
applied to synthesize the mock spectra for stellar rotation, we here adopt the actual 
rotational kernel (eq.[5]) with ui = 0.45, U2 = 0.30, and vsini^ = 4.5 km sec~^. Here 
is the empirical formula after the mock analysis: 



As in the above three cases, we show the OTS formula and the Gaussian formula (37) 
in Figure 5 (d). We adopt a = i;sinz^/1.31 = 3.4 km sec~^, io = 86.7°, and A = 0°. 

The Gaussian formulae in Figure 5 well describe the behavior of simulated results in 
each system even though the actual rotation kernels are not described by simple Gaussians. 
Comparing the three systems except TrES-2 system, in which the stellar velocity is partic- 
ularly small, we notice that the Gaussian formulae for TrES-4 and HD 209458 describes the 
simulation better than that for HD17156. This is presumably due to the fact that apart from 




(45) 




(46) 
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the tiny contribution of the tail, the rotation kernel is better approximated as a Gaussian 
for the star with higher rotational velocity (Appendix F). The fact that the three lines in 
each figure intersect at the same time (close to, but not identical to the central transit time 
time = 0) indicates that while the difference among the OTS formula, the Gaussian formula 
and the simulated result affects the estimate of vsini^, the estimated value of the spin-orbit 
angle is fairly robust as discussed in §2. 

We also notice a fairly big discrepancy between the Gaussian formula and the simulated 
result during the egress of HD17156 (at < time[h] < 1.5) in the upper panel of Figure 5. 
Indeed the OTS formula approximates better the simulated result at the egress phase. In 
order to understand the behavior, we systematically changed the value of io, and found that 
the asymmetry of the velocity anomaly curve between ingress and egress is very sensitive 
to the combination of A and io- This points to a limitation of our Gaussian approximation, 
which might be fixed by using a more precise rotation kernel. 

4. Summary and Discussion 

We have presented an improved analytic formula for the RM effect as measured from 
the cross-correlation method. Our main finding is equation (37) which describes the RM 
velocity anomaly under the Gaussian approximation for the intrinsic line profile and the 
stellar rotation kernel. Unlike the previous approximation (the moment method) adopted by 
OTS, the radial velocity anomaly from the cross-correlation method is explicitly dependent 
both on the thermal and natural broadening width of an individual line profile {(3) and on 
the rotation velocity (a). 

The analytical formula has been compared with the simulated analysis for the Subaru 
HDS routine using the mock spectra constructed from the high-resolution spectrum of the 
Sun. Even though the actual analysis routine attempts to fit simultaneously many different 
lines that are not necessarily approximated by Gaussians, the resulting RM velocity anomaly 
is well described by a form of equation (37). 

The current result explains the previous findings (Winn et al. 2005; Johnson et al. 
2008; Winn et al. 2008) that the OTS formula does not provide a good approximation to 
planetary systems with a relatively high stellar rotation rate. The next step in obtaining a 
more accurate analytic formula for the RM effect would probably be to use a more accurate 
rotation kernel (eq. [5]) instead of the Gaussian approximation (22). For future work we 
plan to pursue this approach, together with the propagation of the resulting precision on the 
misalignment angle A and the rotation velocity wsinz^. 
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In addition to the improvement of the analytic approach, the current result points to a 
possible further refinement of the observational reduction routine for the velocity anomaly 
due to the RM effect. The procedure described in §3.1 is appropriate for the radial veloc- 
ity determination out of transit. During the transit, however, the Doppler-shifted stellar 
spectrum S{\ — 5X) in equation (38) should be replaced with 

>S'star('^ — A) — S'pianct(A — A — A A) 
= 5'star(A — A) — /5'narrowcd template (A — A — AA), (47) 

where S'star(A) is the stellar template spectrum derived by the procedure described in §3.1 
during the out-of-transit data, ^narrowed template (A) is obtained by deconvolving S'star(A) with 
the rotational kernel in principle, and A and AA correspond to the Keplerian plus pecu- 
liar velocity of the star and the radial velocity of the occulted portion of the stellar disk, 
respectively. 

Fitting the observed spectrum in transit with the additional parameters /, A, and AA, 
instead of 5\ alone, one may directly extract the velocity Vp = cAA/Aq after averaged 
over many different lines. Moreover one may reduce the number of fitting parameters by 
estimating / from photometric data during the transit. 

Since the above method is more appropriate to model the spectra of the transiting 
planetary system during transit, one expects that the precision and accuracy of AA, and 
therefore the RM velocity anomaly, may be improved in principle. In reality, however, the 
practical feasibility to implement in the analysis routine depends on the available resolution 
of the spectra as well as their stability. The work toward this improvement is also currently 
under way, which we also hope to present elsewhere in due course. 

The RM effect has become an observationally mature probe of the transiting planetary 
systems. Further improvement in both the accuracy and the precision of the velocity anomaly 
is important to advance the understanding of the formation and evolution processes of such 
systems. We hope that the present analytic approach provides some insights that will enable 
progress in that direction. 

We thank Wako Aoki for useful advice on the spectral line profiles. The radial velocity 
analysis in our simulation was carried out on the computer system at the Astronomy Data 
Center of the National Astronomical Observatory of Japan. N.N. is supported by a Japan 
Society for Promotion of Science (JSPS) Fellowship for Research (PD: 20-8141). Y.S. grate- 
fully acknowledges the support from the Global Scholars Program of Princeton University. 
This work is also supported by JSPS Core-to-Core Program "International Research Net- 
work for Dark Energy". J.N.W. gratefully acknowledges the support of the NASA Origins 
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A. Flux during Complete Transit 



In this Appendix, we compute the flux ratio between the flux coming from the stellar 
disk blocked by the planet and the total flux, and show that equation (16) reproduces the 
OTS formula. 

In the case of quadratic limb darkening (eq. [4]), we can express / in terms of 7, ui, U2, 
and the position of the planet in the stellar disk. The deflnition of / is 



f{x,y) 



I{x',y')dx'dy' 



planet 



(Al) 



I{x', y')dx'dy' 



stellar disk 



where the x- and y-ax.es are deflned so that the y-axis is along the projected rotational axis 
and the origin is at the center of the stellar disk. The integration in the numerator covers 
the portion blocked by planet, while the one in the denominator is performed throughout 
the whole stellar disk. Denoting the radius of the star by Rg, then the cosine in equation (4) 
is written as 



cos^ 



x^ + 



'1- 



(A2) 



Then the denominator of equation (Al) becomes 
-Rs 

dr2TTr 



l-«i<l-^/l-( — 



(A3) 



When the size of the planet is sufficiently small, the numerator of equation (Al) can be 
approximated as 



I{x',y')dx'dy' ~ I{x,y)7rR^,^' 



planet 



2, ,2 



1-ni 1- Vl-P' -^^2 1- Vl-P 



(A4) 
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where p = \J x'^ + y'^/R.s and 7 = Rp/Rg. Substituting equations (A3) and (A4) into equation 
(Al), we obtain 



f{x,y) 



1 - Ml {1 - ^A^} - «2 {1 - vT^y 



1 — ni/3 — U2/Q 
With this expression, equation (16) reduces to 



■7 



Av ~ —Vr. 



1 - ni {1 - vT^} - {1 - ^A^}^ 



(1 - u^/3 - U2/6) - 7^ 



l-Ui 1- Vl-P' - ■«2 1- Vl-P- 



where Vp is expressed by the position of the planet as 



Vp{x,y) 



Qs{x' ,y')x' sini^, dx'dy' 



planet 



dx'dy' 



planet 



(A5) 



(A6) 



(A7) 



The function VLs{x^ y) is the angular velocity of the stellar spin at position (x, y) and i is the 
inclination of the stellar spin axis. For rigidly rotating stars, the planet velocity during a 
complete transit becomes 



Vp = QgX sinij,. 



(A^ 



with Qg being a constant. Equation (A6) with (A8) reproduces the OTS formula for complete 
transit with limb-darkening (eq. [40] in OTS). 



B. During Ingress and Egress 

During an ingress and an egress of transits, the fraction of the portion on the stellar 
disk occulted by the planet changes. The center of the stellar disk and the planet is at the 
origin and at {x,y), respectively. 

Figure 6 shows an example of a configuration during an ingress or an egress. We define 
the radius of the planet as Rp. The angles a and /? are given by 

Rl + r^-Rl Rl + P-Rl 

cosa = — , cos/3 = — , Bl 

ZHgr IRpT 
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Fig. 6. — Schematic illustration of the configuration at an ingress or egress. The center of 
the planet is at {x,y). The dark shaded region is occulted by the planet. 



respectively. The size S of the dark shaded area in Figure 6 (the portion occulted by the 
planet) is 

S = aRl- ]^Rl sin 2a + (3Rl - ^i^J sin 2(3. (B2) 

Thus, the fraction T{r) of the portion on the stellar photosphere occulted by the planet 
becomes 

r(r) = 4^ = - |« - 1 sin2a + - ^ sin2/?)) . (B3) 

During an ingress or an egress, 7^ in equation (A5) must be replaced with r(r) depending 
on the position of the planet. This result is derived in the case of Rs < r, but we can prove 
this also holds in the case of Rg > r. 

The subplanet velocity is no longer expressed by equation (A8). In this case, we must 
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replace equation (A8) with 

/ / ^Igx' sin dx'dy' 

v,{x, y) = -^-^^ , (B4) 

/ / dx'dy' 

where S is the dark shaded region in Figure 6. However, as long as the planet is sufficiently 
small enough, the approximation (A8) provides a good description for subplanet velocity 
(B4). 



C. Derivation of Equation (27) 

We here present the derivation of equation (27) in detail. Below, we implicitly assume 
that the quantities and 7^ represent the broad line originated from the whole photosphere, 
while (5p and 7p indicate the narrow line originated from the portion of the planet. For further 
computation, it is convenient to consider the Fourier transforms of the line proffie functions. 
Note that the Fourier transforms of Gaussian and Lorentzian respectively become 

G'(A) = -^e-^'/^' =^ G{a) = e-^"^"''\ (CI) 
L(A) = --^ =^ L(a)=e-2-^H. (C2) 

Then, the Fourier transform of the Voigt function becomes V{a) = G{cr)L{a). In terms of 
these, equations (25) and (26) are rewritten as 

/oo 
da K(a)e-2™^ - /V;,(a)e-2--(^-^^) , (C3) 
-00 

POO 

^star(A) = - / rfaV;(a)e-2.i.(A)^ ^C4) 



where K(o') and Vp{a) are the Fourier transforms of G(A; /?*)L(A; 7^,) and (^(A; /?p)L(A; 7^), 
respectively. The corresponding Fourier counterpart are 

^transit(a) = -^(a) + /V;(a)e2™^^ (C5) 

^stM = -K(a). (C6) 

Thus, we obtain 

POO 

C((5rm) = / daCia)e^^"''^^, (C7) 
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where 



(C8) 

Now it is easy to incorporate additional effects expressed as a convolution simply by multi- 
plying the Fourier transform of the additional convolution kernel. 

Substituting equations (CI) and (C2) into equation (C8), we obtain 

/oo 
daexp [— 27r^/3^(T^ — 47r7^|a| + 2ma6^^M] 
-oo 

/oo 
daexp [-7r2(/52 + /3^^)a^ - 27r(7, + 7p)|a| + 2mai6nM - AA)] . (C9) 
-oo 

The integral of each term in equation (C9) results in the following integral: 

poo 

^(a,6, c)= / exp{—bx^ + cx)cos{ax)dx. (CIO) 



The derivative of the function ^(a, b, c) with respect to a is 



(a2-c2)/46 



la _ zc)e'"^/2' ( 1 + erfi(^^) ) + c.c. 



2Vb 



(Cll) 



In terms of C,{a, b, c), equation (C9) is expressed as 

C{a) = e (27r5RM, 2n'Pl 4ttj,) - {2n{5^M - AA), Ti\(3l + /^J), 27r(7, + 7^)) . (C12) 

Substituting equation (C12) into equation (17) and using equation (Cll), we finally obtain 
an equation for S-rm'- 



53/2 



- ic)e''"'''^^ ( 1 + erfi( 



,a — ic. 



f 



53/2 



2Vb 

(A-2C)e^^^/2B(l+erfi('^"'^ 



2VB 



where 



a = 27i6km, b = 2n^p^, c = Anj^, 
A = 27r((5RM - AA), 5 = 7r\(3l + (3l), C = 27r(7. + 7^)- 



(C13) 



(CM) 
(C15) 



-27- 



D. Lorentzian profile 

In §2.3, we have derived analytic expression for Gaussian line profiles. In this appendix, 
we consider the opposite simplification that line profile has a Lorentzian form [b = B = 0). 
To derive the radial velocity formula, we use 



6^0 63/2 



iac/2b I ^ „,a — IC, 



(a - ic^e""'''' 1 + erfi( ^ 



(Dl) 



vr [a? + c^Y 



Then, equation (27) or (C13) reduces to 



25rm7* _ J ('^RM - AA)(7^ + 7p) ^^2) 



Employing the similar approximation in §2.3.1, we obtain an approximate expression for 
Lorentzian profile: 

*HM--f^)VAA(l-2(^)1. (D3) 

or 

Again, we obtain the similar expression to equation (37). Note that the factor in front of / 
depends on the widths of Lorentzians, which generally differs from unity. The presence of 
the cubic term in is also a major difference from the result by moment-method. 



E. Velocity Anomaly for the Gaussian Profile Perturbed with Lorentzian 



We discuss a more general case in which the line profile is expressed as Gaussian profile 
with a small contribution of Lorentzian profile. We assume that 7//? is sufficiently small 
(the Gaussian property is dominant) and collect the terms up to the linear order of 7//? in 
equation (27). Noting that 

2z 

eTfi{z) = —= z <^1, 
and defining z = {a — ic)/2Vb, the right hand side of equation (27) becomes 



^Re 






2i V 

















^Re 
b 



z(l - z' 



2i 

1 + —z 

IT 



2a fl 



- 3c' 



'Txh 



(El) 



(E2) 
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Expanding the left hand side of equation (27) in a similar way, and substituting a, b, c and 
A, B, C, we obtain 



2P'^ /AA J 4(7. + 7.) _ AA^ 



or 



This is the formula for the velocity anomaly when absorption lines are Gaussians perturbed 
with Lorentzian profile. Note that equation (E4) reproduces equation (37) when we set 
7* = 7p = 0. 



F. Relation between the Rotation and Gaussian Broadening Kernel by 

Least-squares Fitting 

Here, we derive a scaling factor between the upper-limit wavelength Al of the rotational 
kernel and the dispersion a of the Gaussian broadening kernel by least-squares fitting. Al- 
though replacing a rotational kernel with Gaussian one is a crude approximation, we can 
justify this treatment to some extent (§3.3.2.). 

We define as the residual between the two convolution kernels: 

/oo 
[R{X;\L)-Gi\;P)]''d\ 

/oo /"OO POO 

R\\;\L)d\+ / G\\;p)d\-2 / R{\; \l)G{\; p)d\. (Fl) 
oo J —oo J —oo 

The scaling factor a is define as 

When we fix the rotational velocity (thus A^), a is determined by minimizing of the function 
A^: 

Since the first term in equation (Fl) does not depend on (3, we neglect it in computing the 
derivative in terms of /5. The second term is written as 

oo In/ 

G\X;P)dX = —= = —^. (F4) 
-oo pv27r Alv 2ir 
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The third term in equation (Fl) reduces to 



R{X;\L)G{X;p)d\ 



X, 



+C2 



a 



Cl 



Ci^e-"'/2{/o(aV2) + /i(aV2)} 
+C2-^ |2ae-°' + (-1 + 2a^)^/^ erf(a)| 
+C3^e-"'/2 {«'/o(«V2) + (-1 + «')/i(«V2)} 



Vl — x'^dx 



'(1 — X + C3 



(F5) 



where In{x) is the n-th order modified Bessel function. Substituting equations (6), (7) and 
(8) into equation (F5), we obtain the derivative of A^: 



da 



4 



1 



12(mi + 2M2)e" 



y/2 a^{-6 + 2ui + U2 



a3(-6 + 2Mi+M2) 



{-2a3(-l + Mi+M2)/o(aV2) 



-2a{-2u2 + a^{-l + ui + n2)}/i(aV2) 



v^e"'/2(Mi + 2M2)erf(a^ 



0. 



(F6) 



Substituting the hmb-darkening parameters, we numerically solve the above equation for a. 
We obtain a ~ 1.31 for HD17156 (ui = 0.6, U2 = 0), a ~ 1.30 for TrES-2 (ui = 0.40, 
U2 = 0.30), a ~ 1.32 for TrES-4 {m = 0.46, U2 = 0.31), and a ~ 1.31 for HD209458 
{ui = 0.45, U2 = 0.30). Equation (F6) implies that residuals scales as when we fix 
a. This means that for rapidly rotating star, A^ becomes small and the approximation of 
rotational kernel with Gaussian is validated. 
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